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Abstract. Let V C [0, 1)'^ be a finite point set of cardinality N in an 5- 
dimensional cube, and let / : [0, 1)'^ — > R be an integrable function. A 
QMC integration of / by P is the average of values of / at each point in 
V, which approximates the integration of / over the cube. Assume that V is 
constructed from an F2-vector space P C (^2)^ means of a digital net with 
n-digit precision. As an n-digit discretized version of Josef Dick's method, 
we introduce Walsh figure of merit (WAFOM) WF(P) of P, which satisfies a 
Koksma-Hlawka type inequality, namely, QMC integration error is bounded 
by C's.n I I/I l"WF(P) under n-smoothness of /, where Cs,n is a constant de- 
pending only on S, n. 

We show a Fourier inversion formula for WF(P) which is computable in 
0{nSN) steps. This effectiveness enables us a random search for P with small 
value of WF(P), which would be difficult for other figures of merit such as 
discrepancy. From an analogy to coding theory, we expect that random search 
may find better point sets than mathematical constructions. In fact, a naive 
search finds point sets P with small WF(P). In experiments, we show better 
performance of these point sets in QMC integration than widely used QMC 
rules. We show some experimental evidence on the effectiveness of our point 
sets to even non-smooth integrands appearing in finance. 



1. Introduction 

There is a strong analogy between coding theory and quasi-Monte Carlo (QMC) 
point sets (e.g., see [HI [16] [2] [13]). In coding theory, it seems to be widely 
believed that to find good codes, a random search is easier than mathematical 
explicit constructions. This is supported by the fundamental theorem of Shannon 
showing that a randomly chosen code is good with high probability ^14j. Further 
support is the success of low density parity-check (LDPC) codes [TT]. Thus, it is 
natural to consider a random search for good QMC point sets. An obstruction is 
the lack of a practically computable measure for the "goodness" of QMC point sets. 
In coding theory, the minimum distance is a decisive and computable measure for 
the quality of the codes, so a random search to optimize this value works well, if the 
cardinality of the code words is not large. For QMC point sets, the star-discrepancy 
([12] [^) is hard to compute. Alternative diaphony [18] and dyadic diaphony [8] 
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require 0{SN'^) operations, where S is the dimension and N is the cardinaUty 
of the point sets, which would be practicaUy difficuh for large N for a random 
search. Probably, the most important figure of merit on digital nets is the i- value 
of (t,m, s)-net introduced by Niederreiter ( \V2\ [B]). The i-value gives an upper 
bound on star-discrepancy. Its possible weakness is that t assumes only an integer 
value, and thus may be too coarse to give a tight bound on the star-discrepancy. 

In this paper we propose Walsh figure of merit (WAFOM) as a practically com- 
putable measure of goodness of a QMC point set from a digital net [T^ §4.3]. An 



inversion formula (4.2) computes WAFOM in only 0{nSN) steps, where n is the 
precision (see fj2]) of the point sets. We can easily execute a random search min- 
imizing WAFOM. Our experiment yields QMC point sets which perform better 
than some standard low discrepancy point sets, for QMC integration of a natural 
integrand arising in computational finance. 

Let us briefly recall the notion of QMC integration. Let V C [0, 1)"^ be a point set 
in an S'-dimensional cube with finite cardinality jPI = N, and let / : [0, 1)'^ — >■ M be 
an integrable function. The quasi-Monte Carlo integration by V is an approximation 
value 

1 

fP\ 



(1-1) Mf) E /W 



of the actual integration 



(1.2) /(/) / /(x)dx. 

J[0,1)S 

Thus, the QMC integration error is 

(1.3) Err(/;7') :=!/(/) -/p(/)|. 

A central idea in the theory of QMC is to show an upper bound on this error of 
the form 

Eri{f;V)<V{f)D{V), 

where V{f) is a quantity depending only on / (not V) and D(V) is a quantity 
depending only on V, for a suitable class of integrand functions. Then, one designs 
■p with a small value of D{V), which works for this class of functions. In the 
case of the monumental Koksma-Hlawka inequality, V{f) is the total variation 
of /, D{V) is the star-discrepancy of and the class of functions is those with 
bounded variation. In this case, there are many studies on the construction of 
point sets with small D{V), in particular satisfying a conjectured lower bound of 
0{N-^ (log N)'^-^) (see [H] [Sj for these basics). 

For a given integer a > 1, Dick [3] [4] introduced the notion of a-smooth func- 
tions, the Qf-norm and (he did not give a name so we term it here as) Walsh 
figure-of-merit (WAFOM) WF^^CP), so that 

Err(/;P)<C5,„||/|U-WF„(P) 

holds, using Walsh functions (see [H §4.2]) for a class of point sets (i.e., digital- 
nets). Dick constructed families of point sets with WF^CP) < 0{N-°'{\ogN)°''^). 
A difficulty is that to realize a large value oi a, N becomes too large for practical 
QMC integration. 

Our approach is as follows. In fj2j we fix an integer n called the degree of 
discretization, discretize [0, 1) into 2" intervals, and approximate / by a function 
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/„ which is constant on each smah cube. If we assume Lipschitz continuity of /, the 



approximation error is 0(2~"). In J 3.1 we recall discrete Fourier transformation. 
In { 3.2 we introduce an n-digit discretized version WF(P) of WFq(P), and present 
a result of Dick bounding the integration error by WAFOM in our discretized 
context. In S|4j we give a Fourier inversion formula for WAFOM that reduces the 
computational complexity to 0{nSN). Using this formula, we search for point 
sets with small value of WAFOM by random search. In ^ we show results of 
numerical experiments. We find point sets with small WF(P): the values suggest 
that WF(P) = 0{N^^^^) for some positive constant /3 near one, even in a practical 
range of = 2^*^, 2^^, . . . , 2^^. Experimentation shows that these point sets perform 
better than some standard low discrepancy sequences. 

2. Discretization 

Let us fix a positive integer n, which is called the degree of discretization. To 
simplify the analysis, we discretize / := [0, 1) into a finite set with cardinality 
2". It would be possible to do a similar analysis without such discretization using 
Walsh functions 01 El [S] [H] , but we do not pursue it here and work only on the 
discrete case because we are interested in applications in digital computers where 
such discretization is implicitly done. 

The interval / is decomposed into 2" intervals of equal length [b, 6 + 2^"), where 

(2.1) 6 = 5i2-i + ... + 6„2-", 6i,...,5„ G {0,1}. 

Such a 5 is called an n-digit binary fraction. 

Let F2 — {0, 1} be the two-element field. In F2, arithmetic operations are done 
modulo 2. Let be the space of n-dimensional row vectors: 

F^' = {(5i,...,6„) I 6, eF2}. 

The sum of two such vectors is computed componentwise modulo 2. 

There is a natural identification of a vector (&i, . . . ,6„) with an n-digit binary 



fraction b as in (2.1 ), hence we have an identification of the three sets: Fj , the set 



of n-digit binary fractions, and the set of the 2" intervals, by 

(61, . . . , 5„) = 6i2-i + • • • + 6„2-" ^ h ■■= [b, b + 2'"). 

Let V := {¥2 )^ be the set of S x n matrices with components in F2. As usual, 
B gV is described as B := {bT.j)i<T<s.i<j<n: with bx.j & F2. To B, we associate 

b = (fei, . . . , bs) G with 6t = &t,i2^^ H h6T,n2^", and then an 5-dimensional 

cube Ib ■= lb := X X • • • X Jbg . This gives a discretization of hy V = (F2 )"^. 

Let / : /"^ — K be an integrable function. We define its n-digit discrete approx- 
imation /„ as 

U:V^R, /„(i?) := ^ ^y(x)dx. 

In other words, to B Cz V, we associate an S'-dimensional cube consisting of 
(xi, . . . , xs) such that the first n digits of the binary expansion of xt being equal 
to (6t,i, ■ ■ ■ , ^T,Tt) for each T, and fn{B) is the average value of / over this cube. 
It is easy to see 

(2.2) M E /"(^) ^ / /W'^^- 

' ' Bev 
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In the following, we assume that each point in V has coordinates in n-digit binary 
fractions. Thus, V C corresponds to a subset P C V. 

We define the n-th discretized QMC integration of / by P as 

(2.3) IpAf) ■■=\j^\J2 fniB). 

' ' BeP 

Note that this value is hard to compute in practice, since each fn{B) is an inte- 
gration, but can be approximated by the usual QMC integration /•p(/) as follows. 
Let us define 

||/-/„||oo:= sup \UB)-f{b)\. 
Bev,beiB 

This coincides with the usual supremum norm if /„ is regarded as a (not necessarily 
continuous) function on that takes the (constant) value fn{B) on the hypercube 
Ib- From now on, we use the same notation /„ for this function on . 



Lemma 2.1. The difference between the QMC integral {1.1) and the discretized 



QMC integral (2.3) is bounded by 

\Ivif) - IpAf)\ < ll/-/n||oo. 

In particular, if f is continuous with Lipschitz constant K , namely, if for any 

x,x' e [0,1)^ 

|/(x)-/(x')| <if||x-x'|| 

holds (where the right hand side is the Euclidean norm), then the above value is 
bounded by 

||/-/n|| <ifV^2-". 

Proof. The left hand side of the first inequality is bounded by the average of |/(b) — 
/„(b)| over b G T', hence is bounded by the supremum norm || / — /n||oo- For the 
second inequality, for any b take Ib that contains b. Then, since the diameter of 
Ib is 

|/(b)-/„(b)| < sup \f{^)- f{y)\<K■2-^^^S. 

□ 

By this lemma, when we take n large enough so that K ■ 2^"-\/S' is negligible 
compared to the QMC integration error, then we may identify Ip_n{f) and I-p[f). 
Although the above bound depends on K, in practice, n = 30 would be sufficient for 
typical QMC integration, since QMC integration error would be much larger than 
K2~^^^/S. Further justification for this discretization is that, in single precision 
arithmetic real numbers in / are usually discretized as a 23-digits binary fraction, 
and hence in practical integrations, n > 24 would be sufficiently large. 



After the above discretization, our integration is exactly a finite sum (2.2 1, and 



for a point set P GV, our discretized QMC integration is (2.3). Hence we define 
the n-discretized QMC integration error by the difference: 

(2.4) Err(/;P,n) |/p.„(/) - /(/)|, 



which we are going to use as a proxy for (1.3) with the justification above 
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Remark 2.2. In QMC-intcgration, it would be better to pick the mid point x in 
the hypercube /b and use /(x) as an approximation of the average /n(b) of / over 
/b, instead of using the corner point b. Thus, we should translate V by adding 
(2^"^^, 2^"^^, . . . , 2~"~^). This is reflected in our experiments in ^ 

3. Discrete Fourier transformation 

3.1. Preliminary. This section recalls well-known facts on discrete Fourier trans- 
formations. Recall that V is the set of S* x n matrices with components in F2. A 
subset P C y is an F2-linear subspace if P contains the zero matrix and closed un- 
der summation as F2-matrices (i.e. componentwise sum modulo 2). We henceforth 
assume that P is F2-linear. 

We define a (standard) inner product by 

VxV^¥2, B,A^{B\A):= ^ brjar^j, 

l<T<S,l<j<n 

where the summation in the right hand side is modulo 2. For an F2-linear subspace 
P G V , we define its orthogonal subspace 

P^ ■.= {AeV \ {B\A) = for all B e P}. 

This is a linear subspace of V, and we have 

(3.1) diuiP + dim P^ ^diinV = nS. 

We define 

{B\A) := (-l)(^l-4' e{l,-l}. 

The following lemma is standard (see [15] for a more extensive theoretical treatment 
in the context of coding theory and uniform distribution). 

Lemma 3.1. Let P d V be a linear subspace. Then, 

2-<^l^>-i zfA^P^ 
BeP ^ 

Proof. U A € P-L, then {B\A) = (-l)(^l'4) = 1 for all B e P, which settles the 
first case, li A ^ P-^, then the map P ^ F2 defined by B (P|A) is nonzero 
and F2-linear. This and the dimension formula of linear algebra implies that those 
B G P with {B\A) = has dimension dimP — 1, thus the number of such B is 
|P|/2, and hence is the same with the number of P G P with {B\A) — 1, and the 
sum cancels out. This implies the second case. □ 

For a function 
its discrete Fourier transform / : V — > M is defined by 



Bev 

We have the Fourier expansion formula 

(3.2) /(P)= 5]/(A)(P|A). 

Aev 
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In fact, the right hand side is 



E ^ nB'){B'\A){B\A)) 

Aev \ ' ' B'ev / 



If we consider the summation over A first, the sum of {B'\A){B\A) ~ {B' — B\A) 
is if B' - B ^ y-L = {0}, and \V\ if S' - B = 0. Thus, this sum is /(B). By 
averaging (3.2) over P, we have the following well-known theorem (see [H]). 



Theorem 3.2. (Poisson summation formula) 

Let P C V — (F2)'^ be a linear subspace. For any f : V —?' 
identity 



we have the 



^E/(^) = E/V) r^E(^l^) 



' ' BeP Aev 
We want to know the integration 



BeP 



- E hA). 



AeP^ 



lif) 



^ UB) = /„(0). 



Bev 



Recall that we choose our point set P to be an F2-linear subspace of V . Then, the 



discretized QMC integration ( |2.3[ ) of /„ by P is given by Theorem 3.2 

E /"(^) = E /"(^)- 



(3.3) 

' ' BeP AeP^ 
This implies that the integration error of discretized QMC by P is exactly 

(3.4) r^E/«(^)-/"(0)^ E fr.{A), 

' ' BeP AeP^-{o} 

and an error bound is given by 



(3.5) 



Err(/;P,n) 



E 

AeP^-{o} 



In{A) 



< 



E 

AeP^-{o} 



\In{A)\. 



3.2. Estimation of the Fourier coefficients. Once we have some estimation 
function c{A) > for yl G such that 

\UA)\<Crc{A) 



where C/ is a constant depending only on /, then by (3.5 1 we have an upper 
estimate of the discretized QMC integration error: 



(3.6) 



Err(/;P,n) <C/ 



E 

AeP^-{o} 



c(A). 



Dick gives such a bound for a-smooth functions (see [4] for the definitions of a- 
smooth functions and the norm ||/||a). In our context, the integer a is set to n. 
The following is a very special case of ^^(k) defined in [H §4.1] (the case when the 
base is 2, each coordinate of the integer vector k has a binary expansion of at most 
n digits, and n = a). Note that [21 Propagation rule (2)] (^ Propagation rule (ii)]) 
states that a higher order net which achieves an improved rate of convergence of 
the integration error for a given a achieves an improved rate for all 1 < a' < a. 
Hence it may be possible that an improved rate of convergence can also be achieved 
for a < n. 
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Definition 3.3. For A — {o.T,j)i<T<s,i<j<n G V, define 



l<T<Sa<i<ri 

Note that here each qtj E {0, 1} is considered as an integer, not an clement of F2. 

The following theorem is a part of the results of [4, §4.1]. 

Theorem 3.4. (Josef Dick) Let f be an n-smooth function. There is a constant 
Cs.n depending only on n and S such that 

|/„(^)| <Cs,„||/||„-2-^(^). 

Proof. This is contained in Dick's theorem. We use the following two facts: (1) the 
Walsh coefficient /(k) in |4] coincides with /„(k) when the number of binary digits 
of each component of k do not exceed n, because /„ is the truncation of the Walsh 
series upto the degree n, (2) /„(k) = fn{A) holds for A being the matrix obtained 
as coefficients of binary expansion of k, by definition of Walsh functions. □ 

As in [H §4.2], this theorem and the formula ( |3.5[ ) give the error bound 

Err(/;P,n)<Cs,„||/||„- 2-''^^^- 

AeP^-{a} 

Thus, it is natural to define a kind of "figure of merit" (Walsh figure-of- merit or 
WAFOM) of the point set P: 

Definition 3.5. (WAFOM) 

WF(P) := ^ 2-^(^). 

AGP^-{0} 

Hence we have 
(3.7) Err(/;P,n)<C„^,||/||„.WF(P), 

suggesting that finding P with small values of WF(P) would, a priori, have better 
integration performance. 

Remark 3.6. For several (continuous but not smooth) functions /, we numerically 
compute an approximate value of /„(A) for small n and S, such as n = 5 and 
5 = 3, for each A. For 2"'^' distinct A we observe a tendency of |/„(A)| being 
proportional to 2^^^"*-'. (It is observed for many functions, including the non- 
diffcrcntiable function ( |5.5[ ) defined later. On the other hand, for some special type 
of functions such as the form of /(x) = gi{xi)g2{x2) ■ • ■ gs(xs), the behaviour of 

and is evidence that 



3.4 



|/„(A)| is different.) This supports the bound in Theorem 
the bound ( |3.7[ ) may be adequate in practical applications for some non-smooth 
functions (c.f. |3]), and may suggest a closer relation between the QMC integration 
error and WAFOM. 

So far, we have merely treated a discrete version of Dick's method. Our new 
proposal here is to compute WF(P) by the formula (4.2) in the next section, and 
find P with small WF(P), by some random search of P. 
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For any c : V 
(4.1) 



4. Inversion formula for WF(P) 

Theorem |3.2| (with P and interchanged) gives 



AeP^ BeP 
An exphcit computation below shows the foUowing. 

Theorem 4.1. Let V be the set of S x n matrices with coefficie nts i n ¥2, B 



(br,]) & V , (arj) € V, and let c{A) := 2"^^^) as in Definition 



3.3 Then 



1 



^(^)=M n (i+(-i^'^-2-^). 

^<T<S,l<j<n 

Corollary 4.2. Let P C F be a linear subspace. We have 

^ c{A) = \P^\ J2 KB) = ^ E n (1 + (-1)'"'^2-^-). 



AeP^ BeP 
By subtracting c(0) = 1, we have 



seP i<T<sa<j<n 



(4.2) 



WF(P) 



I PI 



n 



[(l + (_l)bT,,2-J)]_l 



BeP I l<T<S,l<j<n 



This is computable in 0{nSN) steps of arithmetic operations in real numbers, 
where N — \P\. 

Proof. By definition of c, 



^ J2 ciA){B\A) = E 2-2:i<-.,i<.<.^'^-.(P|A) 



1^1 



1^1 



<S,l<j<n i"-T,j (_X)I^l<r<S,l<j<ii '^T.jbTJ 



Aev 



n 

' ' Aev l<T<S,l<j<n 



n 



(l + (-l)''^^^2-J) 



<r<s,i<i<n 



(the last equality translates a sum over 2"'' products into a product of nS of 2 
term sums), which pr oves Theorem 4.1 The first identity in Corollary foll ows from 
|V^| = |P| • |P-L| and (O). Formulal4!2| follows from this and Definition 



3.5 



□ 



A merit of the formula (4.2) is that the number of summation depends only on 
|P|, not |P-L| = 2"S-dimP Definition [3^ 

For QMC, the size |P| = 2'^'"^^ can not exceed a reasonable number of opera- 
tions in a computer, say, dimP < 40, since we need to take average of / over P. 



Thus, (4.2 ) is practically computable, unlike a naive computation of WF(P) from 
Definition 3.5[ which requires an intractable 2"'^~'^'™-^ additions for moderate n 
and 5, say Sn > 80. 
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Remark 4.3. The star-discrepancy of a point set is a standard measure of the 
uniformity of point sets, but it is hard to compute. Thus, only constructions of 
point sets which gives an upper bound on the star-discrepancy are studied, without 



exphcit value, such as (i, to, s)-nets. In contrast, the formula (4.2) makes WAFOM 
computable, which enables a random search. 

As a preceding study, we note that '17] studied optimization of the i-value of 
(i, TO, s)-nets by a randomized (evolutionary) algorithm, and often obtained better 
t-values than the original point sets, such as Sobol point sets and Niederreiter-Xing 
point sets. 

Another type of measures of regularity of a point set is diaphony [TH] (TUl Exer- 
cise 5.27] and dyadic diaphony [8], which are computable in 0{SN'^) steps. 

WAFOM is superior to diaphony in the order of computational complexity, which 
is of practical interest since 0{N'^) steps would be difficult for iterated random 



search for large N. Moreover, WF(P) has a direct error bounding formula (3.7). 

We note that WAFOM can be defined only for digital nets, while star-discrepancy 
and diaphony are measures of uniformity defined on any point sets. 

As stated in the introduction, in coding theory, it is often difficult to construct a 
family of good codes explicitly, and a random search sometimes yield a better result. 
WAFOM and Corollary |4.2| offer a framework for obtaining a good QMC-point set 
by random search in a similar fashion. 

5. Some Numerical Evidence 

5.1. Point set generators. WAFOM is defined for any linear subspace P C V — 
(F2)'^. For fixed n, S and d :— dimP, it is natural to search for P randomly 
by uniform choice of its basis. In this study, we restrict our search to the point 
sets generated by an A/-sequence of a fixed primitive polynomial of degree d, with 
uniform random search of the output transform matrix, which we briefly recall (see 
[7| for basics on Af- sequences). We call this type of point set generator a sequential 



generator, for which WAFOM can be computed in 0{nN) steps (see Proposition 5.1 
below). 

Let (fli, . . . , ad) e Fj. A linear recurring sequence associated to (ai, . . . , a^) is a 
sequence xq, xi, . . . G F2 defined by a recursion in F2: 

(5.1) Xj+d = aiXj+d-i H h adXj, 

with the initial vector {xq, . . . , Xd-i) given. If the polynomial + aif^^^ + ■ ■ • + fld 
is primitive, then the period of the sequence is 2"* — 1 for any nonzero initial vector. 
Such a sequence is called an M-sequence with degree d and exists for any d. The 



set of all the solutions of (5.1) constitutes an F2-vector space of dimension d, since 
the initial vector determines the sequence. Let us fix a nonzero initial vector. 
By the condition on the period, {xi,Xi^i, . . . ,Xi^d-i) is nonzero and distinct for 
< i < 2"^ — 1, and hence assumes all nonzero d-dimensional vectors. Thus , th e 
zero sequence and {xi, Xi+i, . . .) for < z < 2^^ — 1 give all of the solutions to (5.1 1, 



which constitute a d-dimensional F2-vector space. 

For an integer S and < fc < 2'' — 1, define a.n S xd matrix Ck whose (T, j)th en- 
try is Xk+j+T-2, that is: the first row of Ck is (xfc, Xk+i, . . . , Xk+d-i), the second row 
is {xk+i,Xk+2,---,Xk+d), etc., and the last row is {xk+s-i,Xk+s, ■ ■ ■ ,Xk+s+d-2)- 
By the above observation, the set 

W:^{Ck\0<k<2'^-l}U {0} 
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is a d-dimensional sub vector space of the space oi S x d matrices. 

To obtain a d-dimensional subspace P (Z we choose a d x n matrix U of rank 
d, and compute the image WU oi W in V by the muhipUcation of U from the 
right. In other words, we compute S x n matrices CkU for each /c (0 < fc < 2^^ — 1) 
and append the matrix to obtain WU C V. We use WU for our point set P. By 
changing [/, we have various P. 

An advantage of such a construction is in the efficiency for generating points. 
Because the first S~ 1 rows of Ck+iU are the same with the last 5 — 1 rows of CkU, 
to compute the (fc + l)-st point, we need only to compute the last row of Ck+iU . 
To do this, wc compute (xfe+s, Xfc+s+i, • . • ,Xk+s+d-i) by and multiply by U 
from the right. 

Another advantage of sequential generators is in computing WAFOM. In Theo- 
rem |4T] we need to compute 



n 



1<T<S 



W (1 + (-1)^^-2- 

l<j<n 



When br denotes the T-th row of B, this can be written as 



c{B) 



n 

1<T<S 



n '^(^t) 

l<i<n 



where h is defined in an obvious fashion. Since Ck+iU is obtained from CkU 
by removing the first row and attaching the last row, we may obtain c{Ck+iU) 
by simply multiplying c{CkU) by /i(attached row)/ft,(removed row). Hence, if we 
record the values of h for past S rows then only one multiplication and one division 
is necessary. Thus we obtain the following. 

Proposition 5.1. For a point set generated by a sequential generator, WAFOM is 
computable in 0{nN) steps. 

In practice, this multiplication and division accumulate the truncation error. 
One way to avoid this would be to divide the sequence into moderate length sub- 
sequences, and to apply this trick for each subsequence. 

Remark 5.2. We are not sure whether such a choice of point sets harms our ability 
to attain a small value of WAFOM, compared to a random search of P by basis. 

5.2. Finding good point sets. For an application, it is necessary to have a P 
with a small WF(P). As shown in the following experiments, even a naive random 
search turns out to find a good P. 

According to Lemma [2. 1| and the discussion there, we fix ti = 30. We focus on 
the integrand function described in § 5.3 with dimension 5 = 4. Our search for 
a good point set proceeded in two stages for ease. The range of the F2-dimension 
of P is d S {10, 11, ... , 22}. Fix a d in this range. At the first stage, we generate 
a d X d matrix U' of rank d unif ormly at random 5000 times. For each U' , we 
generate WU' C (F2)'^ as in {5.1 then compute WAFOM of WU' with degree of 
discretization d. Thus, we compute 5000 WAFOMS, and identify the best U' . At 
the second stage, we generate an 5x (n — d) matrix uniform at random, concatenate 
it from the right to this best U' to obtain an 5 x n matrix U. Then we compute the 
WAFOM of WU C (F^')'5. We iterate this 2000 times, and take the best P = WU 
with respect to WF(P). 
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Slope of random-search WAFOM: -2.024 




16 
d 



Figure 1. WAFOM (in log scale) for selected point sets. The 
WAFOM of the Sobol sequence is the dashed line above and the 
dotted line below is the WAFOM of random search point sets. The 
thin solid line estimates the rate of convergence. The heavy black 
line has a slope of -1 for reference. 



The coefficients of primitive polynomials (oi, . . . , ad) were generated at random, 
bit by bit, with a bias towards ones: the heuristic being that if many coefficients are 
zero then the point set will satisfy a linear relation with a small number of terms, 
which may be harmful to the effectiveness of QMC. For each d, we used the same 
primitive polynomial of degree d for all matrices U. 

Figure [l] plots the smallest WAFOM values found by this procedure for each d. 
An ordinary least squares estimates of the slope suggests that the best point sets in 
a random search achieve an order of convergence something like 0{N~^~^) where 
/? « 1 in the range of 10 < d < 22. The first 2^^ points of a Sobol sequence are 
a digital net over F2, and therefore form an F2-linear point set, such that we can 
compute the WAFOM; it is included as the dashed line for comparison. The Sobol 
sequence, as well as the Faure and Halton sequences plotted in Figure [2] come from 
the open source C++ library quantlib p]. 

Remark 5.3. Dick [J has shown a construction of a family of point sets with 
WF(P) = 0{N-" (log N)°''^) for arbitrary a > 1. This implies that, if we ex- 
haustively search all of the linear point sets P for every d and plot their smallest 
WAFOM, then we can draw a line with slope —a that is above all the plotted 
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points. The line of smallest WAFOM values in Figure [T] seems to be convex above 
(with an exception at c? = 18), which is in accordance with Dick's result. 

For practical QMC, N is bounded above, and Figure [T] might show the practical 
upper bound of the value a. 

5.3. Performance in QMC integration. To illustrate the effectiveness of this 
method for improving the efficiency of QMC integration, we consider a basic prob- 
lem from computational finance: the pricing of an Asian option. Readers not 



interested in some context can skip directly to (5.5 1, where a problem is formulated 



directly in the form of ( 1.2 1 



The Asian option is a simple example of a path dependent option, meaning that 
the payoff depends upon the path taken to the terminal point. In this case, the 
payoff is some average price of the underlying until maturity. This dependence on 
the value of the underlying at several points in time make its valuation a difficult 
problem, both analytically and computationally. For low dimensional or highly 
specialized cases numerical solution of an associated partial differential equation is 
quite efficient. For extensions beyond this, however, the preferred methodology for 
pricing Asian options is (quasi-) Monte Carlo simulation. 

Here we endeavor to price an option on an asset having a geometric Brownian 
motion dynamics with a current price of Pq, (exponential) volatility tr^, strike K, 
and maturity T year. The riskless discount rate is r. Thus, the price of an Asian 
option with maturity T sampled at times Ti/S, i — 1, . . . , S is 



(5.2) 



where E* denotes expectation with the (logarithmic) drift of P replaced by r and 
x+ = max(a;, 0). Under the risk neutral measure, we can write 



(5.3) 



Pt = Pne 



r-^ ]t+aWt 



where Wt is a standard Wiener process. Inverting the Gaussian distribution, $, we 
can map a quasirandom uniform [0, 1) variate to an increment of Pt by 



(5.4) 



And to a X G [0, 1)"* we can associate a realization of (5.2) using the independent 



increments property of Brownian motion. In particular, for the purposes of QMC 



integration, we write (5.2) as: 



(5.5) / e-^ {\y p^^{:-4-)w,..^s^-\.,) _j\ 



dx. 



As noted in Remark 2.2 our point set is a subset of {0, 1/2", . . . , 1 — 1/2"}"', so 



in our experiment we generate standard normal variates hy x !—> ^^^{x + 1/2 x 2^") 
to center the point within each hypercube. 

The QMC integration errors are presented in Figure [2) where our random search 
method is seen to compare favorably with classical QMC point sets. To evaluate 
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Figure 2. Integration error for the lowest WAFOM generators 
found by the procedure described in Section 5.2 (thin soUd Hne). 
The error of several classical QMC point sets (Halton is dotted, 
Sobol is dashed, Faure is alternating dots and dashes) is included 
for comparison. A heavy black line has slope of -1. 



the "true" option price, we performed QMC integration using a very large classical 
low discrepancy point set. 



6. Conclusion 



We proposed Walsh figure-of-merit (WAFOM, Definition 3.5 1 for those point sets 
P constructed from F2-linear spaces by digital nets, following some earlier work by 
Dick 011]. 

WAFOM is a measure of regularity of the point set, with a mathematical support 



by a Koksma-Hlawka type inequality ( 3.7 ) due to Dick. The WAFOM of a point set 
P C (F2)'^ is quickly computable. Because evaluating this measure requires only 
0{nSN) operations (4.2 1, a random search for generators hav ing good WAFOM 
is a feasible undertaking. Applying even the crude method (j ]5.2| for \P\ = 2'*, 
10 < d < 22, we found point sets with small WAFOM (Figure [T]) with order 
0(7V-i-/3) for ;3 w 1. 

These point sets showed excellent integration performance on an important prob- 
lem for computational finance (Figure |2|, even better than some standard QMC, 
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namely, Sobol, Halton and Faure. In the experiments, the integrand is not differ- 
entiable. This is a supportive evidence of the effectiveness of WAFOM for even non 
smooth functions, which is not covered by Theorem |3.4[ 

Acknowledgements. The authors are thankful to Professor Niederreiter and two 
anonymous referees for valuable comments and for making us aware of some earlier 
literature, and to Professor Hickernell for letting us know [17]. 
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